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The statistical properties of avalanches in a dissipative particulate system under slow shear are 
investigated using molecular dynamics simulations. It is found that the magnitude-frequency dis- 
tribution obeys the Gutenberg-Richter law only in the proximity of a critical density and that the 
exponent is sensitive to the minute changes in density. It is also found that aftershocks occur in this 
system with a decay rate that follows the Modified Omori law. We show that the exponent of the 
magnitude-frequency distribution and the time constant of the Modified Omori law are decreasing 
functions of the shear stress. The dependences of these two parameters on shear stress coincide with 
recent seismological observations [D. Schorlemmer et al. Nature 437, 539 (2005); C. Narteau et al. 
Nature 462, 642 (2009)]. 



Individual earthquakes may be regarded as large scale 
ruptures involving a wide range of structural and compo- 
sitional heterogeneities in the crust. However, the statis- 
tical properties of a population of earthquakes are often 
described by simple power-laws. Among them, two laws 
are ubiquitous and occupy a central position in statistical 
seismology: the Gutenberg-Richter (GR) law |lj and the 
Modified Omori law (MOL) @, 3]. The GR law describes 
the earthquake magnitude-frequency distribution 



P{M) oc 10 



(1) 



where M w is the moment magnitude and b a constant 
with a value around 1 along active fault zones. The MOL 
describes the aftershock occurrence rate 



i(t) oc (t + c)- p , 



(2) 



where t is the elapsed time from the triggering event (the 
so-called mainshock), p a positive non-dimensional con- 
stant with a typical value of 1, and c is a time constant. 
The parameters in these two laws are believed to bear 
some information on the physical state of the crust. In- 
deed, Schorlemmer et al. find that the 6-value of the 
GR law is decreasing going from normal (extension) over 
strike-slip (shear) to thrust (compression) earthquakes 
Narteau et al. also find that the time constant c 
in the MOL has the same dependence on the faulting 
mechanism These two observations indicate that, 
under a simple assumption, b and c are decreasing func- 
tions of shear stress. Although the underlying mechanism 
needs further investigation, this may reflect a common 
time-dependent behavior of fracturing in rocks during the 
propagation of earthquake ruptures and the nucleation of 
aftershocks. 



Because the shear stress along an active fault is not 
directly measurable, a solution to address stress depen- 
dences in earthquake statistics is to analyze models that 
implement a restricted set of physical processes. There 
are indeed many models that resemble seismicity, ranging 
from rock fracture experiments [6-8] to computer simu- 
lations on cellular automata Among them, sheared 
granular media [lol - fhlj fit our purpose perfectly because 
it is a simple representation of a granular fault gouge for 
which both the energy and the stress can be easily de- 
fined. Here, we perform numerical simulations showing 
that, as for real seismicity, avalanches in sheared granu- 
lar matter obeys the GR law and the MOL with b and c- 
values which are decreasing functions of the shear stress. 

Our granular system is made of frictionlcss spheres 
with diameters of d and 0.7c? (the ratio of populations 
is 1 : 1). For the sake of simplicity, we assume that 
the mass M of these particles are the same. We limit 
ourselves to a rather small-size system (N = 1500) for 
computational efficiency. Using the radius and the po- 
sition of particle i, which are denoted by Ri and r^, re- 
spectively, the force between particles i and j is writ- 



ten as F» 



[khij — £n 
, and hij 
Ir 



ij * j 



Ri, 



Here 



% 3 — 1 3 ' 

is the overlap 



length. If Ri + Rj < |rjj|, particles i and j are not in con- 
tact so that the force vanishes. Throughout this study, 
we adopt the units in which d = 1, M = 1, and k = 1. 
We choose C = 2.0, which corresponds to the vanishing 
coefficient of restitution. A constant shear rate 7 is ap- 
plied to the system through the Lees-Edwards boundary 



conditions 15]. Note that under these boundary condi- 



tions the system volume is constant. Thus, the important 
parameters are the shear rate 7 and the packing fraction 
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FIG. 1. Avalanche magnitude-frequency distributions for sev- 
eral densities and shear rates, (a) </> = 0.644. (b) </> — 0.650. 
(c) 7 = 10" 7 except for </> = {0.647, 0.650} (7 = 10" 8 ). 



<f>. A steady state of uniform shear rate 7 can be realized 
starting from a class of special initial conditions. Here, 
we investigate such uniform steady states. 

As observed in many previous works on amorphous 
systems [H-Ei [HHHl , the temporal fluctuations of the 
energy becomes volatile if the shear rate is sufficiently low 
and the density is sufficiently high. From the temporal 
fluctuation of the energy E(t), we define an energy drop 
event as follows. The beginning of an event is defined as 
the time t = t\ at which the energy E(t) starts decreas- 
ing; i. e. , E(ti) = and E{t\) < 0. The definition of the 
end of an event is the opposite; t = £2 at which Efa) = 
and Efa) > 0. Such an energy drop event is referred to 
as an avalanche throughout this paper. We also calcu- 



late the (global) shear stress a(t) using the virial [22 1. 
Important quantities that characterize an avalanche are 
the magnitude M = log 10 [E(ti) — E(t 2 )} + 11, the initial 
stress a = <r(ti), the stress drop Act = <j{ti) — ct(£ 2 ), and 
the duration T = t<i — t\. Note that, with this defini- 
tion of magnitude, the GR law reads P(M) oc i(r 2 / 3bM , 
as M w ~ 2/3\og 10 [E(t 1 ) — Efa)] + const. Hereafter, we 
use (3 = 2/36 instead of b. Note also that here the spa- 
tial information of avalanches is discarded. In this case, 
one might overlook simultaneous avalanches occurring at 
different places. However, this is unlikely because the 
present system is small (i. e. the characteristic length is 
approximately of 9d). 

First, we discuss the nature of the avalanche 
magnitude-frequency distribution. Figures [T] (a) and 
(b) show these distributions at several shear rates for 
4> = 0.644 and (f> — 0.650, respectively. The distribution 
is independent of the shear rate at sufficiently low shear 
rates. The characteristic shear rate below which the dis- 
tribution is rate-independent depends on the density and 
may be interpreted as the inverse of the structural re- 
laxation time. In the density range investigated here, 
the characteristic shear rate is a decreasing function of 
the density. For example, it is approximately 10~ 6 for 
4> = 0.644 and 10 -8 for <p = 0.650. Hereafter, we discuss 
such rate-independent behaviors by choosing sufficiently 



low shear rates. 

Figure []] (c) shows rate-independent magnitude- 
frequency distributions at several densities. They may 
be regarded as the GR law in the proximity of a crit- 
ical density ((f) ~ 0.644), whereas they no longer obey 
the GR law at higher densities. At much lower densi- 
ties {(f> < 0.642), an avalanche is no longer observed. It 
should be noted that the (3- value seen in the proximity of 
a critical density is sensitive to the minute changes in den- 
sity. The /3-value is a decreasing function of the density 
ranging from 0.47 to 0.78. We obtain the smallest value 
(0.47) at the largest density {(/> = 0.645) and the largest 
value (0.78) at the smallest density {(f) = 0.642). The 
(3- value at intermediate densities {<fr = {0.643, 0.644}) is 
approximately 0.64. This range of values is comparable 
to that of earthquakes for which the /3-value varies from 
0.50 along normal faults (i. e. , extensional regime, low 
stress) to 0.73 along reverse faults (i. e. , compressional 
regime, high stress). 

We remark that distribution functions of other quanti- 
ties also obey power laws. Among them, we find that the 
distribution functions of the duration T exhibit power- 
law tails, the exponent of which is approximately 3.0 ir- 
respective of the density. We also find that the exponent 
for the stress drop distribution is twice larger than that 
for energy drop. 

As similar power law behaviors have been observed in 
some granular systems, it is interesting to compare the 
present result with other studies. A major interest is the 
/3-value: 0.36 to 0.95 [11], 0.82 to 0.89 [14], and 1 [H have 
been reported, whereas the GR law is not observed in an 
experiment [l2j . Although such non- universality may be 
perplexing to statistical physicists, it may be rather ra- 
tional in view of the present study; namely, the GR law 
holds only in the proximity of a critical density and the 
/3-value depends on both the density and the stress level. 
Thus, the /3-value may not be an analogue for critical 
exponents. As an evidence of such non-universality, one 
can further illustrate a wide range of /3-values obtained 
in various amorphous systems (l6l - [2l1 |. Thus, there may 
exist many unknown ingredients that control the /3-value, 
and further investigation is required for the unified un- 
derstanding of all these exponents in general amorphous 
systems. 

To discuss the effect of shear stress on the /3-value, we 
introduce the shear stress ct at the beginning of an event 
as an additional argument to the magnitude-frequency 
distribution. Then, P(M, ct) is the conditional proba- 
bility of observing a magnitude M avalanche under the 
(global) shear stress value ct. For convenience, ct is in- 
tegrated in a certain interval Si = [10~ 4 / 2 , 10~ l / 2+0 - 5 ] 
(i = {1, 2, • • • , 9}). Thus, we obtain the following dis- 
tribution function: P t (M) = J s daP(M,a). Figured] 
shows the behaviors of Pi(M) at cf> = 0.643. We can see 
that the probability of observing a larger avalanche in- 
creases as the shear stress increases. More importantly, 
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FIG. 2. Avalanche magnitude-frequency distribution for dif- 
ferent ranges of shear stress values (see text for the definition). 
The /3-value decreases at higher shear stress. The volume frac- 
tion is <p — 0.643 and the shear rate is 7 = 10~ 7 . 



the distribution function at each stress level develops the 
tail, with the /3-value being a decreasing function of the 
shear stress. We confirm that this stress dependence is in- 
dependent of the density in the range 0.642 < 4> < 0.645. 
In addition, the shear stress dependence of the /3-value 
is qualitatively the same as that in rock fracture experi- 
ments pj. 

Because aftershocks result from changes of stress 
induced by a mainshock, we disregard the smaller 
avalanches for which M < and consider ranges 
of magnitude for mainshocks [-M™- n , M^f ax ] and after- 
shocks [M^ in , M^ ax ]. In practice, a magnitude M G 
[M^f in , M^f ax ] event occuring at time tyi is not a main- 
shock if there is at least one avalanche of the same 
or higher magnitude range in the time interval [tM — 
AT, tM + AT]. Thus, we only consider isolated main- 
shocks and, taking sufficiently large AT-values, we also 
avoid overlapping aftershock sequences that belong to 
different mainshocks. Then, all M < M^ in avalanches 
that follow a magnitude M M mainshock within the time 
interval [tn^M + AT] are regarded as aftershocks. Fi- 
nally, each aftershock is characterized by its magnitude 
M A and the elapsed time r since the mainshock. 

In order to reduce artifacts related to event detectabil- 
ity, we disregard larger magnitude range for mainshocks 
and smaller magnitude range for aftershocks. In ad- 
dition, the magnitude range of mainshocks are chosen 
within the intermediate magnitude range in which the 
GR law holds. Using this strict methodology, we con- 
siderably reduce the number of events in our artificial 
catalogues. Therefore, in all ranges of magnitudes, af- 
tershocks are stacked with respect to the time of their 
mainshocks to finally end-up with a single aftershock se- 
quence. 

Figure [3] (a) shows the aftershock frequencies for sev- 
eral choices of [M^ in , M^J. The magnitude range 
of aftershocks is chosen as [1,3]. We find that the 
MOL, Eq. ([2]), holds with the exponent p ~ 1 irre- 
spective of the magnitude range of mainshocks. Fig- 
ures [3] (b) and (c) show the density dependence of 



FIG. 3. Aftershock frequencies as functions of the elapsed 
time from the mainshock. The shear rate is 7 = 10" 7 . Af- 
tershocks are stacked for several (> 10) mainshocks and the 
number of aftershocks for each data set is typically larger than 
1000. (a) = 0.644, M A G [1,3], AT = 2 x 10 4 . The dashed 
line is proportional to 1/t. (b) M M G [3,4], M A G [1,3], 
AT = 5 x 10 3 . (c) M M e [4, 5], M A € [2, 4], AT = 1 x 10 4 . 
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FIG. 4. Aftershock frequencies at (/> = 0.644 for several stress 
ranges; a — [exp(i — 10), exp(i — 9)] (i = 0, 1, 2, 3). The magni- 
tude ranges are M M G [3, 4], and M A G [1, 3]. Each frequency 
is normalized as n(0) = 1 so that the stress dependence of c 
is apparent. The inset shows the stress dependence of the 
c-value. Circles: = 0.644, M M G [3,4], and M A G [1,3]. 
Squares: = 0.645, M M G [3, 4], and M A G [1, 3]. Diamonds: 
(t> = 0.644, M M G [2, 3], and M A G [1, 2]. 



aftershock frequencies. We test several densities and 
shear rate (</> = {0.642, 0.643, 0.644, 0.645} and 7 = 
{10~ 6 , 10~ 7 , 10~ 8 }) to find that the exponent p ~ 1 is 
insensitive to the values of these parameters. This ro- 
bustness makes a contrast to the sensitivity of the expo- 
nent P in the GR law to the density. We also find that 
the time constant c is independent of the density. More 
importantly, we find that the time constant c is indepen- 
dent of the shear rate. This indicates that aftershocks 
are not driven by the additional shear strain applied af- 
ter a mainshock; rather, they are caused by the intrinsic 
relaxation dynamics. Thus, we do not expect aftershocks 
in the quasi-static shear deformation, in which a system 
fully relaxes to a stable configuration after an avalanche. 

Next, we show that the time constant c depends on the 
shear stress. We define for aftershocks the stress range 
[c m in, Cmax] an d select only aftershocks that belong to 
this stress range and the magnitude range [M^ in , M^ ax ]. 



4 



The aftershock frequencies are shown in Figure |H in 
which the MOL (with p = 1) holds clearly and c is 
a decreasing function of shear stress. We estimate the 
c- values by fitting the data with Aj(r + c) using the 
maximum-likelihood method. As shown in the inset of 
Figure HI the c- value has a negative dependence on the 
shear stress. We confirm this stress dependence for two 
densities {(p = {0.644, 0.645}) and for two magnitude 
ranges (M M G [3,4], M A G [1,3] and M M G [4,5], 
M A G [2,4]). This negative shear-stress dependence of 
the c-value for aftershocks is similar to those suggested 
by seismological observations and various models of af- 
tershock production. As in models based on friction, 
damage mechanics, or static fatigue, our numerical ex- 
periment shows that a power-law relaxation regime may 
be delayed for an amount of time which is inversely pro- 
portional to the global shear stress. 

To conclude, two fundamental laws in earthquake 
statistics (GR and MOL) are also relevant to avalanches 
in sheared granular matter. Most importantly, the stress 
dependences of the parameters in these laws are the same 
as those observed in real seismicity. This coincidence 
may be rather remarkable in view of the simpleness of 
the simulation model and the complexity of the crust. 
A quantitative result in the present study, that the time 
constant c in MOL is inversely proportional to the ap- 
plied shear stress, may provide seismologists with a new 
tool for inferring the stress level in the crust. It may also 
be relevant to prevention of silo avalanches, because one 
can infer the stress level in silos by observing a time se- 
ries of acoustic emissions. To address such applications, 
the present result must be tested in more realistic and 
complex systems that are relevant to seismogenic zones 
and various industrial situations. 

The author (TH) acknowledges helpful discussions 
with H. Matsukawa and M. Otsuki. 
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